# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)# tablelibrary(gt)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar anglelibrary(oce)# data wranglinglibrary(tidyr)library(dplyr)library(purrr)library(forcats)library(lubridate)
2 Setting up custom function
2.1windrose
Code
windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet = F) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, 7, 20), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin lableslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet) { data2$origin <-unique(data_to_plot$origin) }# set up max value maxvalue <-35# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="white",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet) { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }
3 Import data
Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.
Code
# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")
Code
# get solar angle when data locationsolar_angle_inter <- data_dive %>%filter(!is.na(Lat)) %>%mutate(solar_angle =sunAngle(date, Long, Lat)$altitude) %>%select(id, DiveNumber, solar_angle)# merge that with data_divedata_dive <-merge(data_dive, solar_angle_inter,by =c("id", "DiveNumber"),all.x = T)# add day-night# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/data_dive <- data_dive %>%mutate(day_night =if_else(solar_angle <-18, "night", "day"))
Let’s add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, i.e. if this difference is above xxx m then is not considered as a benthic dive.
Code
data_dive <- data_dive %>%# change DiveType to include condition on difference# between max depth and bathymetry (<50m)mutate(DiveType_50 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >50, 0, DiveType),DiveTypeName_50 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >50,"Transit", DiveTypeName ),# same with <100mDiveType_100 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >100, 0, DiveType),DiveTypeName_100 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >100,"Transit", DiveTypeName ),# same with <150mDiveType_150 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >150, 0, DiveType),DiveTypeName_150 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >150,"Transit", DiveTypeName ),# same with <200mDiveType_200 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >200, 0, DiveType),DiveTypeName_200 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >200,"Transit", DiveTypeName ) )
4 Maps
Code
# check if background_ggoceanmaps existif (file.exists("../export/background_ggoceanmap.rds")) { trip <-readRDS("../export/background_ggoceanmap.rds")} else {# using ggOceanMaps trip <-basemap(limits =c(170, -110, 20, 59),bathymetry =TRUE,shapefiles ="Arctic",rotate =TRUE,grid.col =NA )# Make the graticules: lims <-attributes(trip)$limits graticule <- sf::st_graticule(c(lims[1], lims[3], lims[2], lims[4]),crs =attributes(trip)$proj,lon =seq(-180, 180, 45),lat =seq(-90, 90, 10) )# Plot trip <- trip +geom_sf(data = graticule, color ="grey50") trip$layers <- trip$layers[c(1, 3, 2)]# save resultsaveRDS(trip, "../export/background_ggoceanmap.rds")}
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="All Benthic Dives",subtitle =paste(nrow(df_kernel_dens), "dives") )
Figure 1: Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_50 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<50",subtitle =paste(nrow(df_kernel_dens), "dives") )
Figure 2: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_100 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<100",subtitle =paste(nrow(df_kernel_dens), "dives") )
Figure 3: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_150 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<150",subtitle =paste(nrow(df_kernel_dens), "dives") )
Figure 4: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
Code
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_200 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<200",subtitle =paste(nrow(df_kernel_dens), "dives") )
Figure 5: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of dives per individual per hour" )
Figure 11: Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time
Code
# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
statistic
p.value
method
alternative
62027.5
8.80423e-09
Wilcoxon rank sum test with continuity correction
two.sided
Code
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of dives per individual per hour" )
Figure 12: Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time
Code
# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
statistic
p.value
method
alternative
94857.5
9.608421e-07
Wilcoxon rank sum test with continuity correction
two.sided
Code
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of dives per individual per hour" )
Figure 13: Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time
Code
# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
statistic
p.value
method
alternative
94412.5
1.916629e-06
Wilcoxon rank sum test with continuity correction
two.sided
Code
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of dives per individual per hour" )
Figure 14: Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time
Code
# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
statistic
p.value
method
alternative
94853.5
9.670154e-07
Wilcoxon rank sum test with continuity correction
two.sided
Code
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date) ) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day",y ="Average number of dives per individual per hour" )
Figure 15: Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time
Code
# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")
Mann Whitney / Wilcoxon rank sum test on benthic dives rate
statistic
p.value
method
alternative
94946.5
1.710967e-06
Wilcoxon rank sum test with continuity correction
two.sided
7 Boxplot of the number of benthic dive (day vs. night)
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 1000)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 16: Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 17: Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 18: Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 19: Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time
data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim =c(0, 250)) +labs(x ="Time of the day",y ="Number of benthic dives per individual" )
Figure 20: Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time
---title: "Benthic Dives Investigations"author: - name: "Astarte Brown" - name: "Conner Hale" - name: "Joffrey Joumaa"date: "`r invisible(Sys.setlocale(locale = 'C')); format(Sys.Date(), format = '%B %d, %Y')`"format: html: toc: true toc-location: left number-sections: true smooth-scroll: true code-fold: true code-tools: true df-print: paged fig-align: "center" highlight-style: arrow self-contained: trueexecute: echo: true cache: false warning: falsetheme: light: flatly dark: darklyknitr: opts_chunk: message: false rownames.print: false tidy: styler---# Import library```{r}# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)# tablelibrary(gt)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar anglelibrary(oce)# data wranglinglibrary(tidyr)library(dplyr)library(purrr)library(forcats)library(lubridate)```# Setting up custom function## `windrose````{r}windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet = F) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, 7, 20), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin lableslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet) { data2$origin <-unique(data_to_plot$origin) }# set up max value maxvalue <-35# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="white",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet) { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }```# Import dataLet's load `data_dive`, *i.e.* the output of `data_wrangling.qmd`, and filter only on animals leaving from Ano Nuevo.```{r}# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")``````{r}# get solar angle when data locationsolar_angle_inter <- data_dive %>%filter(!is.na(Lat)) %>%mutate(solar_angle =sunAngle(date,Long,Lat)$altitude) %>%select(id, DiveNumber, solar_angle)# merge that with data_divedata_dive =merge(data_dive, solar_angle_inter,by =c("id", "DiveNumber"), all.x = T)# add day-night# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/data_dive <- data_dive %>%mutate(day_night =if_else(solar_angle <-18, "night", "day"))```Let's add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, *i.e.* if this difference is above xxx m then is not considered as a benthic dive.```{r}data_dive = data_dive %>%# change DiveType to include condition on difference# between max depth and bathymetry (<50m)mutate(DiveType_50 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >50, 0, DiveType),DiveTypeName_50 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >50,"Transit", DiveTypeName ),# same with <100mDiveType_100 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >100, 0, DiveType),DiveTypeName_100 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >100,"Transit", DiveTypeName ),# same with <150mDiveType_150 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >150, 0, DiveType),DiveTypeName_150 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >150,"Transit", DiveTypeName ),# same with <200mDiveType_200 = data.table::fifelse(DiveType ==3& (abs(bathy) - Maxdepth) >200, 0, DiveType),DiveTypeName_200 = data.table::fifelse( DiveTypeName =="Benthic"& (abs(bathy) - Maxdepth) >200,"Transit", DiveTypeName ) )```# Maps```{r}# check if background_ggoceanmaps existif (file.exists("../export/background_ggoceanmap.rds")) { trip <-readRDS("../export/background_ggoceanmap.rds")} else {# using ggOceanMaps trip <-basemap(limits =c(170, -110, 20, 59),bathymetry =TRUE,shapefiles ="Arctic",rotate =TRUE,grid.col =NA )# Make the graticules: lims <-attributes(trip)$limits graticule <- sf::st_graticule(c(lims[1], lims[3], lims[2], lims[4]),crs =attributes(trip)$proj,lon =seq(-180, 180, 45),lat =seq(-90, 90, 10) )# Plot trip = trip +geom_sf(data = graticule, color ="grey50") trip$layers <- trip$layers[c(1,3,2)]# save resultsaveRDS(trip, "../export/background_ggoceanmap.rds")}```:::{.panel-tabset}## All```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-all}#| fig-cap: "Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="All Benthic Dives",subtitle =paste(nrow(df_kernel_dens), "dives"))```## diff<50```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_50 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-50}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<50",subtitle =paste(nrow(df_kernel_dens), "dives"))```## diff<100```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_100 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-100}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<100",subtitle =paste(nrow(df_kernel_dens), "dives"))```## diff<150```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_150 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-150}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<150",subtitle =paste(nrow(df_kernel_dens), "dives"))```## diff<200```{r}# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat) & DiveTypeName_200 =="Benthic") %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName)``````{r fig-map-200}#| fig-cap: "Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park."#| fig-height: 5#| fig-asp: 0.6# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plottrip +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# titlelabs(title ="Benthic Dives with bathy-max_depth<200",subtitle =paste(nrow(df_kernel_dens), "dives"))```:::# Windrose:::{.panel-tabset}## All```{r fig-windrose-all}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(3000, 6000, 9000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff<50```{r fig-windrose-50}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<50) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_50 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff<100```{r fig-windrose-100}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<100) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_100 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff<150```{r fig-windrose-150}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<150) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_150 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```## Diff<200```{r fig-windrose-200}#| fig-cap: "Circular histogram plots displaying the time when they perform benthic dives (diff<200) across their foraging trip."#| # let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic")# figurewindrose_benthic <- data_windrose %>%group_by(id) %>%filter(DiveTypeName =="Benthic") %>%mutate(time =as.numeric(format(date_tz, format ="%H"))) %>%group_by(time) %>%summarise(nb_ind_hour =n()) %>%mutate(origin ="Benthic") %>%windrose(., grid =c(1000, 2000, 3000), legend_position ="top", facet = T)# displaywindrose_benthic``````{r}# perform testtest_benthic =r.test(conversion.circular(circular( data_windrose %>%group_by(id) %>%filter(DiveTypeName_200 =="Benthic") %>%mutate(time =hour(date_tz) +minute(date_tz)/60+second(date_tz)/(60*60)) %>%pull(time),units ="hours" ),units ="radians"))# convert back into human reading hourconv_test =conversion.circular(circular(test_benthic$r.bar, units ="rad"), units ="hours")# printprint(paste0("The average time is ", round(conv_test[[1]],1), " (p-value: ", round(test_benthic$p.value, 1), ")"))```:::# Boxplot of the benthic dive rate (day vs. night):::{.panel-tabset}## All```{r fig-boxplot-all}#| fig-cap: "Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time"#| data_plot = data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of dives per individual per hour")``````{r}# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")```## diff<50```{r fig-boxplot-50}#| fig-cap: "Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of dives per individual per hour")``````{r}# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")```## diff<100```{r fig-boxplot-100}#| fig-cap: "Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of dives per individual per hour")``````{r}# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")```## diff<150```{r fig-boxplot-150}#| fig-cap: "Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of dives per individual per hour")``````{r}# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")```## diff<200```{r fig-boxplot-200}#| fig-cap: "Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%mutate(hour = lubridate::hour(date),day =as.Date(date)) %>%group_by(id, day, hour, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop") %>%group_by(id, day_night) %>%summarise(nb_average_benthic_dive =mean(nb_benthic_dive), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +geom_boxplot(outlier.shape =NA) +labs(x ="Time of the day", y ="Average number of dives per individual per hour")``````{r}# Mann Whitney / Wilcoxon rank sum test on benthic dive ratebroom::tidy(wilcox.test( data_plot %>%filter(day_night =="day") %>%pull(nb_average_benthic_dive), data_plot %>%filter(day_night =="night") %>%pull(nb_average_benthic_dive))) %>%gt() %>%tab_header(title ="Mann Whitney / Wilcoxon rank sum test on benthic dives rate")```:::# Boxplot of the number of benthic dive (day vs. night):::{.panel-tabset}## All```{r fig-boxplotbis-all}#| fig-cap: "Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,1000))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff<50```{r fig-boxplotbis-50}#| fig-cap: "Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_50 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff<100```{r fig-boxplotbis-100}#| fig-cap: "Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_100 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff<150```{r fig-boxplotbis-150}#| fig-cap: "Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_150 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```## diff<200```{r fig-boxplotbis-200}#| fig-cap: "Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time"#| data_plot <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# keep all benthic divesfilter(DiveTypeName_200 =="Benthic") %>%group_by(id, day_night) %>%summarise(nb_benthic_dive =n(), .groups ="drop")data_plot %>%ggplot(aes(x = day_night, y = nb_benthic_dive)) +geom_boxplot(outlier.shape =NA) +coord_cartesian(ylim=c(0,250))+labs(x ="Time of the day", y ="Number of benthic dives per individual")``````{r}broom::tidy(glm( nb_benthic_dive ~ day_night,family ="quasipoisson",data = data_plot)) %>%gt() %>%fmt_number(columns =2:5,decimals =2,suffixing =TRUE) %>%tab_header(title ="nb_benthic_dive ~ day_night")```:::